INF-495, v0.01, Claudio Torres, ctorres@inf.utfsm.cl. DI-UTFSM

Textbook: Computational Mathematical Modeling, An Integrated Approach Across Scales by Daniela Calvetti and Erkki Somersalo.

Chapter 2



In [13]:

    
import numpy as np
import scipy.sparse.linalg as sp
import sympy as sym
from scipy.linalg import toeplitz

import ipywidgets as widgets
from ipywidgets import IntSlider

import matplotlib.pyplot as plt
%matplotlib inline
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter

plt.style.use('ggplot')

import matplotlib as mpl
mpl.rcParams['font.size'] = 14
mpl.rcParams['axes.labelsize'] = 20
mpl.rcParams['xtick.labelsize'] = 14
mpl.rcParams['ytick.labelsize'] = 14
sym.init_printing()

from scipy.integrate import odeint
from ipywidgets import interact

2.1 Simple compartment models

Example 1:



In [21]:

    
# Example 1, numerically.
def f_example_1(t,y,c,k1,k2):
    X = y[0]
    Y = y[1]
    
    Xp = c(t)-k1*X
    Yp = k1*X-k2*Y
    
    return np.array([Xp,Yp])

# Initial condition
y0 = np.array([500, 600])
a  = 1200
t1 = 12*60
t2 = 19*60
tau = 30
c  = lambda t: a*(np.exp(-((np.mod(t,24*60)-t1)**2)/(2*tau**2)) \
                +np.exp(-((np.mod(t,24*60)-t2)**2)/(2*tau**2))) \
                /(np.sqrt(2*np.pi*tau**2))

T_digest = 80
T_metabolic = 240
k1 = 1/T_digest
k2 = 1/T_metabolic

T = 2*24*60

# time where we want your solution
t = np.linspace(0, T, 100)
sol = odeint(f_example_1, y0, t, args=(c,k1,k2), tfirst=True)

plt.figure()
plt.plot(t, sol[:,0], 'b', label='X(t): Energy in gastric channel')
plt.plot(t, sol[:,1], 'r', label='Y(t): Energy in blood')
plt.legend(loc='right', bbox_to_anchor=(2, 0.5))
plt.xlabel('t')
plt.ylabel('kcal')
plt.grid(True)
plt.show()



In [18]:

    
plt.figure()
T = 5*24*60
t = np.linspace(0, T, 1000)
plt.plot(t,c(t),'b')
plt.grid(True)
plt.plot()









    Out[18]:





$\displaystyle \left[ \right]$

2.1.1.- Simple population models

Adding competition for resources: Logistic growth



In [27]:

    
# Example 1, numerically.
def f_spm_logistic_growth(t,y,r,K):
    X = y
    Xp = r*(1-X/K)*X
    return Xp

# Initial condition
r = 1
K = 1000

T = 10

# time where we want your solution
t = np.linspace(0, T, 100)

plt.figure()
for y0 in np.linspace(0,2000,20):
    sol = odeint(f_spm_logistic_growth, y0, t, args=(r,K), tfirst=True)
    plt.plot(t, sol, 'b')
plt.xlabel('t')
plt.ylabel('population size')
plt.grid(True)
plt.show()

Logistic growth with harvesting



In [41]:

    
def f_logistic_growth_harvesting(t,y,r,K,h):
    X = y
    Xp = r*(1-X/K)*X-h
    return Xp

# Initial condition
r = 1
K = 1000
h = 100
print(K*r/4,h)

T = 8

Xplus = K/2*(1+np.sqrt(1-4*h/(K*r)))
Xminus = K/2*(1-np.sqrt(1-4*h/(K*r)))
print(Xplus,Xminus)

# time where we want your solution
t = np.linspace(0, T, 100)

plt.figure()
N=20
y0s = np.zeros(N+1)
y0s[0] = Xminus-0.5
y0s[1:] = np.linspace(Xminus,1200,N)
for y0 in y0s:
    sol = odeint(f_logistic_growth_harvesting, y0, t, args=(r,K,h), tfirst=True)
    plt.plot(t, sol, 'b')

plt.plot(t, sol*0+Xplus, 'r')
plt.plot(t, sol*0+Xminus, 'r')
plt.plot(t, sol*0, 'r--')
plt.xlabel('t')
plt.ylabel('population size')
plt.grid(True)
plt.show()









    



250.0 100
887.2983346207417 112.7016653792583



In [42]:

    
plt.figure()
f = lambda X: -(r/K)*X**2+r*X-h
Xs = np.linspace(0,Xplus+200,100)
plt.plot(Xs, f(Xs), 'b')
plt.plot(Xs, f(Xs)*0, 'r')
plt.axvline(x=Xminus,color='k',label='$X_-$')
plt.axvline(x=Xplus,color='g',label='$X_+$')
plt.xlabel('X')
plt.legend(loc='right', bbox_to_anchor=(1.3, 0.5))
plt.grid(True)
plt.show()

Models with delay

Example 2

Differential equations with delay are not currently supported by Scipy, fortunately there are the following options:

For Python: pydelay, http://pydelay.sourceforge.net
For R: "PBSddesolve: Solver for delay differential equations", https://github.com/pbs-software/pbs-ddesolve
For Python using PBSddesolve: https://github.com/hensing/PyDDE
(*) For Python: https://github.com/Zulko/ddeint



In [43]:

    
# See: https://github.com/Zulko/ddeint
from ddeint import ddeint

r = 2         # Growth rate
K = 100       # Carrying capacity
tau = 1       # Lag of crowding
X_hist = 50   # Constant history
T = 30        # Time simulation

def values_before_zero(t):
    return X_hist

def f_delay1(Y, t):
    return r*(1-Y(t - tau)/K)*Y(t)

t = np.linspace(0, T, 1000)
sol = ddeint(f_delay1, values_before_zero, t)

fig, ax = plt.subplots(1, figsize=(4, 4))
ax.plot(t, sol)
plt.show()



In [45]:

    
r = 1.4       # Growth rate
K = 150       # Carrying capacity
tau1 = 1      # Lag of birth
tau2 = 5      # Lag of crowding
X_hist = 100  # Constant history
T = 300       # Time simulation
beta = 0.5    # Rate of death

def values_before_zero(t):
    return X_hist

def f_delay2(Y, t):
    return r*Y(t-tau1)-beta*Y(t)-(r/K)*Y(t - tau2)*Y(t)

t = np.linspace(0, T, 2000)
sol = ddeint(f_delay2, values_before_zero, t)

fig, ax = plt.subplots(1, figsize=(4, 4))
ax.plot(t, sol)
plt.show()

2.2.- Interacting population models

2.2.1.- Predator-prey models

Example 3



In [48]:

    
# Example 1, numerically.
def f_example_3(t,y,eta,gamma,eps):
    X = y[0]
    Y = y[1]
    
    Xp = eta*X*(1-Y/(eta/gamma))
    Yp = -eps*Y*(1-X/(eps/delta))
    
    return np.array([Xp,Yp])

# Initial condition

eta = 8-1/2
eps = 1/5
gamma = 1/100
delta = eps*gamma/(20*eta)

X_eq = eps/delta
Y_eq = eta/gamma

y0 = np.array([0.9*X_eq, 1.2*Y_eq])

T = 15

# time where we want your solution
t = np.linspace(0, T, 100)
sol = odeint(f_example_3, y0, t, args=(eta,gamma,eps), tfirst=True)

X_output = sol[:,0]
Y_output = sol[:,1]

fig, ax1 = plt.subplots()
color = 'tab:blue'
ax1.set_xlabel('time [years]')
ax1.set_ylabel('Prey population', color=color)
ax1.plot(t, X_output, color=color)
ax1.tick_params(axis='y', labelcolor=color)
ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis
color = 'tab:green'
ax2.set_ylabel('Predator popolation', color=color)  # we already handled the x-label with ax1
ax2.plot(t, Y_output, color=color)
ax2.tick_params(axis='y', labelcolor=color)
fig.tight_layout()  # otherwise the right y-label is slightly clipped
plt.show()

plt.figure()
plt.plot(t,X_output,'b')
plt.plot(t,Y_output,'r')
plt.grid(True)
plt.show()



In [9]:

    
plt.figure()
plt.plot(X_output,Y_output,'k')
plt.xlabel('prey')
plt.ylabel('predator')
plt.plot(X_eq,Y_eq,'b.',markersize=20)
plt.plot(y0[0],y0[1],'r.',markersize=20)
plt.text(X_eq,Y_eq,'  equilibrium state')
plt.text(y0[0],y0[1],'  initial state')
plt.show()



In [10]:

    
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(X_output, Y_output, t,'b')
ax.view_init(elev=40,azim=230)
plt.xlabel('prey')
plt.ylabel('predator')
ax.set_zlabel('t')
plt.show()

2.2.2.- SIR models in epidemiology



In [52]:

    
# Example 1, numerically.
def f_SIR(t,y,alpha,beta):
    S = y[0]
    I = y[1]
    R = y[2]
    
    Sp = -alpha*I*S
    Ip = alpha*I*S-beta*I
    Rp = beta*I
    
    return np.array([Sp,Ip,Rp])

# Initial condition
N = 1e6
alpha = 1e-6
beta = 1/3

S0 = 9e5
I0 = 1e5
y0 = np.array([S0, I0, N-S0-I0])

T = 20

# time where we want your solution
t = np.linspace(0, T, 100)
sol = odeint(f_SIR, y0, t, args=(alpha,beta), tfirst=True)

plt.figure()
plt.plot(t, sol[:,0], 'b', label='S(t)')
plt.plot(t, sol[:,1], 'r', label='I(t)')
plt.plot(t, sol[:,2], 'g', label='R(t)')
plt.legend(loc='right', bbox_to_anchor=(1.3, 0.5))
plt.xlabel('t')
plt.ylabel('number of individuals')
plt.grid(True)
plt.show()

2.2.3.- Endemic diseases



In [53]:

    
# Example 1, numerically.
def f_SIR_endemic(t,y,alpha,beta,sigma,N):
    S = y[0]
    I = y[1]
    R = y[2]
    
    Sp = -alpha*I*S+sigma*N-sigma*S
    Ip = alpha*I*S-beta*I-sigma*I
    Rp = beta*I-sigma*R
    
    return np.array([Sp,Ip,Rp])

# Initial condition
N = 1e6
alpha = 1e-6
beta = 1/3
sigma = 1/50

S0 = 9e5
I0 = 1e5
y0 = np.array([S0, I0, N-S0-I0])

T = 150

# time where we want your solution
t = np.linspace(0, T, 100)
sol = odeint(f_SIR_endemic, y0, t, args=(alpha,beta,sigma,N), tfirst=True)

plt.figure()
plt.plot(t, sol[:,0], 'b', label='S(t)')
plt.plot(t, sol[:,1], 'r', label='I(t)')
plt.plot(t, sol[:,2], 'g', label='R(t)')
plt.legend(loc='right', bbox_to_anchor=(2, 0.5))
plt.xlabel('t')
plt.ylabel('number of individuals')
plt.grid(True)
plt.show()

Project 1

https://github.com/projectmesa/mesa

https://towardsdatascience.com/introduction-to-mesa-agent-based-modeling-in-python-bcb0596e1c9a



In [ ]: